clear all
close all
%testImage=imread('C:\Users\Michael Ferry\Desktop\Bridget\ymc003\t000308xy2c1.tif');
%testImage=imread('C:\Users\Michael Ferry\Desktop\Bridget\ymc003\t000001xy2c1.tif');
testImage=imread('C:\Users\Michael Ferry\Desktop\WillFilter\t000304xy3c1_shape.tif');
outimage='C:\Users\Michael Ferry\Desktop\WillFilter\t000304xy3c1_shape_Test.tif';

figure, imagesc(testImage)
title('Original image')
colormap('gray')

%seBall = strel('ball',3,3);


PSF = fspecial('gaussian',5,5);
Blurred = imfilter(testImage,PSF,'symmetric','conv');
%PSF = fspecial('disk',2);

% seBall = strel('ball',3,3);
% testImageErode=imdilate(Blurred,seBall);

%Create a simulated blur in the image and add noise.
%BlurredPre = imfilter(testImage,PSF,'symmetric','conv');

%Find edges
%PSF_H = fspecial('sobel');
%PSF_V = PSF_H';
%BlurredH = double(imfilter(testImage,PSF_H,'symmetric','conv'));
%BlurredV = double(imfilter(testImage,PSF_V,'symmetric','conv'));
%Blurred = sqrt(BlurredH.^2 + BlurredV.^2);

%imnoise(Blurred,'gaussian',0,V);
figure,imagesc(Blurred);
%figure,imagesc(testImageErode);
title('blurred, eroded image')
%image 2


colormap('gray')
seDisk = strel('disk',3);
testImageBS=imtophat(testImage,seDisk);
%testImageBS=imtophat(Blurred,seDisk);

figure,imagesc(testImageBS);
title('Tophat filtered');

%testImageBS=imtophat(testImage,seDisk);

%testImageAdaptedPre = adapthisteq(testImageBS,'Range','full','Distribution','exponential');
testImageAdapted=histeq(testImageBS);

%imshow(testImageAdapted);
%testImageAdapted = medfilt2(testImageAdaptedPre);
level=graythresh(testImageAdapted);
bw = im2bw(testImageAdapted,level);
bw2 = medfilt2(bw);

figure,imshow(bw2)
title('Image after thresholding and median filter')

seDiskClose = strel('disk',6);
bw3Close=imclose(bw2,seDiskClose);
figure,imshow(bw3Close)
title('Closed image')


blobsizeLow=100;
blobsizeHigh=2000;

bw3 = bwareaopen(bw2,blobsizeLow)+bwareaopen(~bw3Close,blobsizeHigh);

figure,imshow(bw3)
title('Image after blob removal')
imwrite(bw3,outimage,'jpg');

seDiskClose = strel('disk',1);
bw3Close=imclose(bw3,seDiskClose);
figure,imshow(bw3Close)
title('Closed image')

bw4 = bwdist(bw3);

 figure, imshow(bw4,[],'InitialMagnification','fit')
 title('Distance transform of ~bw3')

bw4 = -bw4;
bw4(bw3) = -Inf;




bw5=~watershed(bw4);
%bw5=~watershed(bw3,4);
%seBall = strel('ball',3,3);
%bw4=imdilate(bw3,seBall);
%hline = strel('line',3,0);
%vline = strel('line',3,90);
%sqSe = strel('square',3);
%bw4=imdilate(bw3,hline);
%bw5=imdilate(bw4,vline);

%bw5=imdilate(bw3,sqSe);

figure,imshow(bw5)
title('Watershed of filtered image')
%image5
% rgb = label2rgb(bw5,'jet',[.5 .5 .5]);
% figure, imshow(rgb,'InitialMagnification','fit')
% title('Watershed transform of D')

bw6=logical((bw5==1)+bw3);

figure,imshow(bw6)
title('original image plus watershed lines')

blobsizeLow2=25;
blobsizeHigh2=1500;
bw6_blobR = ~bwareaopen(~bw6,blobsizeLow2)+bwareaopen(~bw6,blobsizeHigh2);
figure,imshow(bw6_blobR)
title('blob removal')

% CurrentObjects = regionprops(bwconncomp(~bw6,4),'Area','Centroid','Eccentricity','Orientation','MajorAxisLength','MinorAxisLength','PixelIdxList');
% %CurrentObjects = regionprops(bwconncomp(~bw3Close,4),'Area','Centroid','Eccentricity','Orientation','MajorAxisLength','MinorAxisLength','PixelIdxList');
% dim=[1040,1392];
% bw7= true(dim(1),dim(2));
% areaCutoffLow=50;
% areaCutoffHigh=2000;
% eccCutoff= 0.9;
% 
% for j=1:length(CurrentObjects)
%     
%     object=j;
%     if(CurrentObjects(object).Area > areaCutoffLow && CurrentObjects(object).Area < areaCutoffHigh)
%         if(CurrentObjects(object).Eccentricity < eccCutoff)
%             center_1 = CurrentObjects(object).Centroid;
%             %o = regionprops (L, 'orientation');
%             orientations = cat(1, CurrentObjects(object).Orientation);
%             phi_1 = orientations(1)* (pi/180)  ; %convert to radians
%             %phi_1 = - phi_1;
%             
%             %a = regionprops (L, 'majorAxisLength');
%             majorAxisLengths = cat(1, CurrentObjects(object).MajorAxisLength);
%             a_1 = majorAxisLengths(1)/2 ;
%             
%             %b = regionprops (L, 'minorAxisLength');
%             minorAxisLengths = cat(1, CurrentObjects(object).MinorAxisLength);
%             b_1 = minorAxisLengths(1)/2 ;
%             
%             t= (0 : 0.01 : 2*pi);
%             
%             % Draw Ellipse
%             x_1 = round(center_1(1)+a_1*cos(phi_1)*cos(t)-b_1*sin(phi_1)*sin(t));
%             y_1 = round(center_1(2)+a_1*sin(phi_1)*cos(t)+b_1*cos(phi_1)*sin(t));
%             %idx= y_1*dim(1)+ x_1;
%             
%             for i=1:length(x_1)
%                 if(round(x_1(i)) > 0 && round(y_1(i)) > 0)
%                 %if(any(idx < 0))
%                     %testOut(round(x_1(i)),round(y_1(i)))=0;
%                     bw7(round(y_1(i)),round(x_1(i)))=0;
%                 %else
%                     %testOut(idx)=0;
%                 end
%             end
%         end
%     end
% end
% bw8=~imfill(~bw7,'holes');
% figure
% imshow(bw7)
% figure
% imshow(bw8)

% center=CurrentObjects(20).Centroid;
% phi=CurrentObjects(20).Orientation;
% a=CurrentObjects(20).MajorAxisLength;
% b=CurrentObjects(20).MinorAxisLength;
% 
% t=[0:0.1:2*pi()];
% 
% 
% x = center(2)+a*cos(phi)*cos(t)-b*sin(phi)*sin(t); 
% y = center(1)+a*sin(phi)*cos(t)+b*cos(phi)*sin(t); 



% [M N]=size(testImageBS);
% phaseRESHAPE=reshape(double(testImageBS),M*N,1);
% figure
% hist(phaseRESHAPE,25)

%%bw3 = bwareaopen(bw2,blobsize);
        
%phaseRESHAPE_sort=sort(phaseRESHAPE);   